Chapter 10 Changes after FMT compared to acclimation

<<<<<<< HEAD
load("data/data_27022025.Rdata")
load("data/calculations_28022025.Rdata")

10.1 What is the trend of the microbiota in each type after FMT?

10.1.1 Alpha diversity

label_map <- c(
  "Control" = "Cold-control", 
  "Treatment" = "Cold-intervention",
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity"
)

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point %in% c("FMT1","Acclimation", "FMT2")) %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
#  facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
      scale_fill_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA50","#d57d2c50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=10))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2") 

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(5.7695)  ( log )
Formula: richness ~ type * time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   757.7    783.7   -367.8    735.7       68 
=======
load("data/data_27022025.Rdata")
load("data/calculations_28022025.Rdata")

10.1 What is the effect of FMT and temperature on the GM after 1 week compared to acclimation?

10.1.1 CI vs CC

10.1.1.1 Alpha diversity

label_map <- c(
  "Control" = "Cold-control", 
  "Treatment" = "Cold-intervention",
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity"
)
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(type %in% c("Control","Treatment") & time_point %in% c("FMT1","Acclimation")) %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
#  facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=10))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>%
  filter(type!="Hot_control") 

Richness

Modelq0GLMMNB <- MASS::glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)

Call:
MASS::glm.nb(formula = richness ~ type * time_point, data = alpha_div_meta, 
    trace = TRUE, init.theta = 2.532718573, link = log)

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    3.9512     0.2275  17.368   <2e-16 ***
typeTreatment                 -0.1372     0.3132  -0.438    0.661    
time_pointFMT1                -0.3225     0.3140  -1.027    0.304    
typeTreatment:time_pointFMT1   0.4500     0.4435   1.015    0.310    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(2.5327) family taken to be 1)

    Null deviance: 37.769  on 33  degrees of freedom
Residual deviance: 36.407  on 30  degrees of freedom
AIC: 327.2

Number of Fisher Scoring iterations: 1

              Theta:  2.533 
          Std. Err.:  0.630 

 2 x log-likelihood:  -317.205 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
 type      emmean    SE  df asymp.LCL asymp.UCL
 Control     3.79 0.157 Inf      3.48      4.10
 Treatment   3.88 0.157 Inf      3.57      4.18

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast            estimate    SE  df z.ratio p.value
 Control - Treatment  -0.0878 0.222 Inf  -0.396  0.6921

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 

Neutral

Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta) 
summary(Modelq1n)

Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)

Residuals:
    Min      1Q  Median      3Q     Max 
-20.699  -6.754   1.842   6.297  15.828 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    21.373      3.716   5.751  2.8e-06 ***
typeTreatment                  -3.974      5.107  -0.778    0.443    
time_pointFMT1                 -4.014      5.107  -0.786    0.438    
typeTreatment:time_pointFMT1   11.478      7.223   1.589    0.123    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.51 on 30 degrees of freedom
Multiple R-squared:  0.08999,   Adjusted R-squared:  -0.001011 
F-statistic: 0.9889 on 3 and 30 DF,  p-value: 0.4113

Phylogenetic

Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta) 
summary(Model_phylo)

Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.01194 -0.49654  0.06198  0.55787  1.90952 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.02586  0.1608  
Number of obs: 79, groups:  individual, 27

Fixed effects:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                      3.9089     0.1696  23.046   <2e-16 ***
typeHot_control                  0.5568     0.2283   2.439   0.0147 *  
typeTreatment                   -0.1233     0.2297  -0.537   0.5916    
time_pointFMT1                  -0.3140     0.2167  -1.449   0.1473    
time_pointFMT2                   0.1487     0.2173   0.684   0.4937    
typeHot_control:time_pointFMT1  -0.1155     0.2983  -0.387   0.6987    
typeTreatment:time_pointFMT1     0.4309     0.3068   1.404   0.1602    
typeHot_control:time_pointFMT2  -0.3876     0.2988  -1.297   0.1945    
typeTreatment:time_pointFMT2     0.4093     0.2990   1.369   0.1710    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typHt_cntrl -0.741                                                         
typeTretmnt -0.721  0.535                                                  
tim_pntFMT1 -0.677  0.503  0.497                                           
tim_pntFMT2 -0.702  0.520  0.507  0.532                                    
typH_:_FMT1  0.496 -0.671 -0.362 -0.727 -0.389                             
typTr:_FMT1  0.482 -0.357 -0.664 -0.707 -0.378  0.515                      
typH_:_FMT2  0.515 -0.684 -0.370 -0.388 -0.730  0.516     0.276            
typTr:_FMT2  0.498 -0.370 -0.686 -0.384 -0.719  0.280     0.511    0.524   
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
 type        emmean     SE  df asymp.LCL asymp.UCL
 Control       3.85 0.1050 Inf      3.65      4.06
 Hot_control   4.24 0.0996 Inf      4.05      4.44
 Treatment     4.01 0.1030 Inf      3.81      4.21

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE  df z.ratio p.value
 Control - Hot_control     -0.389 0.144 Inf  -2.711  0.0184
 Control - Treatment       -0.157 0.145 Inf  -1.083  0.5246
 Hot_control - Treatment    0.232 0.143 Inf   1.627  0.2342

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
 time_point  emmean     SE  df asymp.LCL asymp.UCL
 Acclimation   4.05 0.0938 Inf      3.87      4.24
 FMT1          3.84 0.0949 Inf      3.66      4.03
 FMT2          4.21 0.0895 Inf      4.03      4.38

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE  df z.ratio p.value
 Acclimation - FMT1    0.209 0.123 Inf   1.700  0.2050
 Acclimation - FMT2   -0.156 0.121 Inf  -1.286  0.4028
 FMT1 - FMT2          -0.365 0.122 Inf  -2.997  0.0077

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

Neutral

Modelq1n <- nlme::lme(neutral ~ type*time_point, 
               random = ~ 1 | individual,
               data = alpha_div_meta)
summary(Modelq1n)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  563.2664 587.9998 -270.6332

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    4.551076 9.160157

Fixed effects:  neutral ~ type * time_point 
                                    Value Std.Error DF   t-value p-value
(Intercept)                     21.343055  3.603126 46  5.923483  0.0000
typeHot_control                 23.258307  4.960549 24  4.688656  0.0001
typeTreatment                   -3.944072  4.960549 24 -0.795088  0.4344
time_pointFMT1                  -3.983722  4.472618 46 -0.890691  0.3777
time_pointFMT2                   2.374489  4.472618 46  0.530895  0.5980
typeHot_control:time_pointFMT1 -14.877819  6.216964 46 -2.393101  0.0208
typeTreatment:time_pointFMT1    11.253000  6.325237 46  1.779064  0.0818
typeHot_control:time_pointFMT2 -14.775497  6.216964 46 -2.376642  0.0217
typeTreatment:time_pointFMT2    20.108022  6.216964 46  3.234380  0.0023
 Correlation: 
                               (Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typeHot_control                -0.726                                                         
typeTreatment                  -0.726  0.528                                                  
time_pointFMT1                 -0.663  0.481  0.481                                           
time_pointFMT2                 -0.663  0.481  0.481  0.534                                    
typeHot_control:time_pointFMT1  0.477 -0.649 -0.346 -0.719 -0.384                             
typeTreatment:time_pointFMT1    0.469 -0.340 -0.638 -0.707 -0.378  0.509                      
typeHot_control:time_pointFMT2  0.477 -0.649 -0.346 -0.384 -0.719  0.518     0.272            
typeTreatment:time_pointFMT2    0.477 -0.346 -0.649 -0.384 -0.719  0.276     0.509    0.518   

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.1324340 -0.7088533  0.0456110  0.6569569  1.9681295 

Number of Observations: 79
Number of Groups: 27 
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
 type        emmean   SE df lower.CL upper.CL
 Control       20.8 2.36 26     16.0     25.7
 Hot_control   34.2 2.33 24     29.4     39.0
 Treatment     27.3 2.36 24     22.4     32.2

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                estimate   SE df t.ratio p.value
 Control - Hot_control     -13.37 3.31 24  -4.038  0.0013
 Control - Treatment        -6.51 3.33 24  -1.952  0.1461
 Hot_control - Treatment     6.86 3.31 24   2.073  0.1170

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation   27.8 2.01 24     23.6     31.9
 FMT1          22.6 2.01 24     18.4     26.7
 FMT2          31.9 1.97 24     27.9     36.0

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Acclimation - FMT1     5.19 2.55 46   2.034  0.1156
 Acclimation - FMT2    -4.15 2.52 46  -1.646  0.2372
 FMT1 - FMT2           -9.34 2.52 46  -3.703  0.0016

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

Phylogenetic

Model_phylo <- nlme::lme(phylogenetic ~ type*time_point, 
               random = ~ 1 | individual,
               data = alpha_div_meta)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  269.7976 294.5311 -123.8988

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.4986226 1.145917

Fixed effects:  phylogenetic ~ type * time_point 
                                   Value Std.Error DF   t-value p-value
(Intercept)                     5.256251 0.4407611 46 11.925398  0.0000
typeHot_control                 1.266332 0.6064636 24  2.088060  0.0476
typeTreatment                   0.283689 0.6064636 24  0.467776  0.6442
time_pointFMT1                 -0.837008 0.5590601 46 -1.497169  0.1412
time_pointFMT2                 -0.481515 0.5590601 46 -0.861293  0.3935
typeHot_control:time_pointFMT1 -1.410952 0.7774021 46 -1.814958  0.0761
typeTreatment:time_pointFMT1   -0.545877 0.7906304 46 -0.690432  0.4934
typeHot_control:time_pointFMT2 -0.585100 0.7774021 46 -0.752635  0.4555
typeTreatment:time_pointFMT2    0.056099 0.7774021 46  0.072162  0.9428
 Correlation: 
                               (Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typeHot_control                -0.727                                                         
typeTreatment                  -0.727  0.528                                                  
time_pointFMT1                 -0.676  0.492  0.492                                           
time_pointFMT2                 -0.676  0.492  0.492  0.533                                    
typeHot_control:time_pointFMT1  0.486 -0.663 -0.353 -0.719 -0.383                             
typeTreatment:time_pointFMT1    0.478 -0.348 -0.652 -0.707 -0.377  0.509                      
typeHot_control:time_pointFMT2  0.486 -0.663 -0.353 -0.383 -0.719  0.517     0.271            
typeTreatment:time_pointFMT2    0.486 -0.353 -0.663 -0.383 -0.719  0.276     0.509    0.517   

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.81722823 -0.52018811 -0.06481484  0.55835817  2.23628404 

Number of Observations: 79
Number of Groups: 27 
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
 type        emmean    SE df lower.CL upper.CL
 Control       4.82 0.280 26     4.24     5.39
 Hot_control   5.42 0.276 24     4.85     5.99
 Treatment     4.94 0.280 24     4.36     5.52

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE df t.ratio p.value
 Control - Hot_control     -0.601 0.393 24  -1.527  0.2963
 Control - Treatment       -0.120 0.396 24  -0.304  0.9505
 Hot_control - Treatment    0.481 0.393 24   1.221  0.4525

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE df lower.CL upper.CL
 Acclimation   5.77 0.245 24     5.27     6.28
 FMT1          4.28 0.245 24     3.78     4.79
 FMT2          5.12 0.241 24     4.62     5.61

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE df t.ratio p.value
 Acclimation - FMT1    1.489 0.319 46   4.666  0.0001
 Acclimation - FMT2    0.658 0.316 46   2.085  0.1042
 FMT1 - FMT2          -0.831 0.316 46  -2.635  0.0302

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

10.1.1.1 Beta diversity

Number of samples used

<<<<<<< HEAD
samples_to_keep_post7 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2") %>% 
  select(Tube_code) %>% 
  pull()
subset_meta_post7 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2")

subset_meta_post7$time_point<-as.factor(subset_meta_post7$time_point)
subset_meta_post7$type<-as.factor(subset_meta_post7$type)
length(samples_to_keep_post7)
[1] 79

Richness

richness_post7 <- as.matrix(beta_q0n$S)
richness_post7 <- as.dist(richness_post7[rownames(richness_post7) %in% samples_to_keep_post7,
               colnames(richness_post7) %in% samples_to_keep_post7])
betadisper(richness_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
=======
samples_to_keep_post3 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post3 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control")
subset_meta_post3$type_time<-interaction(subset_meta_post3$type, subset_meta_post3$time_point)
length(samples_to_keep_post3)
[1] 34

Richness

richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
               colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     2 0.06454 0.032268 4.9753    999   0.01 **
Residuals 76 0.49290 0.006485                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
<<<<<<< HEAD
              Control Hot_control Treatment
Control                 0.0160000     0.697
Hot_control 0.0119513                 0.002
Treatment   0.6817089   0.0035106          
adonis2(richness_post7 ~ type*time_point,
        data = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))),
        permutations = 999,
        strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
======= Control Treatment Control 0.484 Treatment 0.49714
adonis2(richness_post3 ~ type*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_b9npyawa3sc60nzpd9sg, .table th.tinytable_css_b9npyawa3sc60nzpd9sg { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_375734g7wmtqf6o2enbc, .table th.tinytable_css_375734g7wmtqf6o2enbc { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 2 3.120806 0.11961754 5.546155 0.001
time_point 2 1.703776 0.06530413 3.027874 0.001
type:time_point 4 1.570884 0.06021051 1.395852 0.001
Residual 70 19.694403 0.75486782 NA NA
Total 78 26.089870 1.00000000 NA NA
<<<<<<< HEAD
subset_meta_post7_arrange <- column_to_rownames(subset_meta_post7, "Tube_code")
subset_meta_post7_arrange<-subset_meta_post7_arrange[labels(richness_post7),]

pairwise <- pairwise.adonis(richness_post7, subset_meta_post7_arrange$type, perm=999)
pairwise%>%
  tt()
=======
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(richness_post3),]

pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_tbs2p2aupt3txp1lbmlr, .table th.tinytable_css_tbs2p2aupt3txp1lbmlr { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_a6pij8s035i9eixvhixm, .table th.tinytable_css_a6pij8s035i9eixvhixm { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 10.7820564 2.401526 0.04582931 0.001 0.003 *0.3657243 1.123239 0.06966584 0.248 1.000
Control vs Hot_control 1 2.6071063 9.057490 0.15081366 0.001 0.003 *
Treatment vs Hot_control 1 1.2773606 4.350043 0.07859150 0.001 0.003 *
pairwise <- pairwise.adonis(richness_post7, subset_meta_post7_arrange$time_point, perm=999)
pairwise%>%
  tt()
<<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 1 1.1985353 3.576695 0.06675842 0.001 0.003 *
Acclimation vs FMT2 1 0.7814316 2.468536 0.04616802 0.003 0.009 *
FMT1 vs FMT2 10.5620083 1.802651 0.03413940 0.011 0.033 .0.5615418 1.729004 0.10335366 0.017 0.102

Neutral

<<<<<<< HEAD
neutral_post7 <- as.matrix(beta_q1n$S)
neutral_post7 <- as.dist(neutral_post7[rownames(neutral_post7) %in% samples_to_keep_post7,
               colnames(neutral_post7) %in% samples_to_keep_post7])
betadisper(neutral_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
=======
neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
               colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
<<<<<<< HEAD
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.06738 0.033689 3.8309    999  0.024 *
Residuals 76 0.66833 0.008794                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.2060000     0.173
Hot_control 0.1985046                 0.005
Treatment   0.1611083   0.0069718          
adonis2(neutral_post7 ~ type*time_point,
        data = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))),
        permutations = 999,
        strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))) %>% pull(individual),
        by="terms") %>%        
        broom::tidy() %>%
        tt()
======= Df Sum Sq Mean Sq F N.Perm Pr(>F) Groups 1 0.019686 0.0196859 2.8068 999 0.11 Residuals 32 0.224435 0.0070136 Pairwise comparisons: (Observed p-value below diagonal, permuted p-value above diagonal) Control Treatment Control 0.119 Treatment 0.10361
adonis2(neutral_post3 ~ type,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual)) %>%        
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_z8wsgq48zbrihjwlnh7p, .table th.tinytable_css_z8wsgq48zbrihjwlnh7p { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_dj3n5nfx6r2y9emr2fzn, .table th.tinytable_css_dj3n5nfx6r2y9emr2fzn { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 2 3.705559 0.14953527 7.609355 0.001
time_point 2 2.217656 0.08949199 4.553951 0.001
type:time_point 4 1.813192 0.07317011 1.861693 0.001Model 1 0.4361022 0.04050485 1.350872 1
Residual 70 17.044094 0.68780262 NA NA
Total 78 24.780501 1.00000000 NA NA
<<<<<<< HEAD
subset_meta_post7_arrange <- column_to_rownames(subset_meta_post7, "Tube_code")
subset_meta_post7_arrange<-subset_meta_post7_arrange[labels(neutral_post7),]
pairwise <- pairwise.adonis(neutral_post7, subset_meta_post7$type, perm=999)
pairwise%>%
  tt()
=======
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(neutral_post3),]
pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_iagsoasjtf3ipwhzulgy, .table th.tinytable_css_iagsoasjtf3ipwhzulgy { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_a18ysd90ugpwhxe6heg0, .table th.tinytable_css_a18ysd90ugpwhxe6heg0 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 11.027801 3.463893 0.064789390.2759858 0.9928976 0.06208366 0.483 1.000
Control.Acclimation vs Control.FMT1 1 0.8153894 3.0970603 0.171136100.001 0.003 *
Control vs Hot_control 13.120025 12.118158 0.19199163 0.001 0.003 *1.1809241 4.4856470 0.24265567 0.002 0.012 .
Treatment vs Hot_control 1 1.394947 5.015972 0.08954539 0.001 0.003 *
Treatment.Acclimation vs Treatment.FMT1 1 1.3127773 4.6298256 0.23585668 0.001 0.006 *
Control.FMT1 vs Treatment.FMT1 1 0.6051778 2.2508491 0.13047758 0.015 0.090
<<<<<<< HEAD
pairwise <- pairwise.adonis(neutral_post7, subset_meta_post7_arrange$time_point, perm=999)
pairwise%>%
  tt()
=======

Phylogenetic

phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
               colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.04697 0.046974 2.8076    999  0.101
Residuals 32 0.53540 0.016731                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Control Treatment
Control               0.096
Treatment 0.10357          
adonis2(phylo_post3 ~ type,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual)) %>%        
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_khvi88bgx62bf4888wka, .table th.tinytable_css_khvi88bgx62bf4888wka { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_5y8xbfi1thjnoay3akf7, .table th.tinytable_css_5y8xbfi1thjnoay3akf7 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 11.6914939 5.559526 0.10006431 0.001 0.003 *0.03710337 0.01720754 0.5602822 1
Acclimation vs FMT2 1 0.9608387 3.212301 0.05925409 0.001 0.003 *
FMT1 vs FMT2 1 0.6282036 2.174336 0.04089070 0.010 0.030 .
<<<<<<< HEAD

Phylogenetic

phylo_post7 <- as.matrix(beta_q1p$S)
phylo_post7 <- as.dist(phylo_post7[rownames(phylo_post7) %in% samples_to_keep_post7,
               colnames(phylo_post7) %in% samples_to_keep_post7])
betadisper(phylo_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
=======

NMDS

beta_richness_nmds_post3 <- richness_post3 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post3, by = join_by(sample == Tube_code)) %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Treatment"),
                       labels=c("Cold-Cold", "Cold-Hot"),
                       values=c("#4477AA","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_neutral_nmds_post3 <- neutral_post3 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post3, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Treatment"),
                       labels=c("Cold-Cold", "Cold-Hot"),
                       values=c("#4477AA","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_phylogenetic_nmds_post3 <- phylo_post3 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post3, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Treatment"),
                       labels=c("Cold-Cold", "Cold-Hot"),
                       values=c("#4477AA","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")
ggarrange(beta_richness_nmds_post3, beta_neutral_nmds_post3, beta_phylogenetic_nmds_post3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

10.1.2 Functional differences

CC from acclimation to FMT1

element_gift %>%
  filter(time_point %in% c("FMT1","Acclimation") & type == "Control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 11 × 6
   trait  p_value p_adjust Domain       Function                                        Element         
   <chr>    <dbl>    <dbl> <chr>        <chr>                                           <chr>           
 1 B0204 0.00370    0.0478 Biosynthesis Amino acid biosynthesis_Serine                  Serine          
 2 B0205 0.000165   0.0117 Biosynthesis Amino acid biosynthesis_Threonine               Threonine       
 3 B0703 0.00370    0.0478 Biosynthesis Vitamin biosynthesis_Niacin (B3)                Niacin (B3)     
 4 B0805 0.00156    0.0370 Biosynthesis Aromatic compound biosynthesis_Indole-3-acetate Indole-3-acetate
 5 D0102 0.000329   0.0156 Degradation  Lipid degradation_Fatty acid                    Fatty acid      
 6 D0212 0.000987   0.0350 Degradation  Polysaccharide degradation_Arabinan             Arabinan        
 7 D0304 0.00247    0.0438 Degradation  Sugar degradation_D-Arabinose                   D-Arabinose     
 8 D0306 0.00156    0.0370 Degradation  Sugar degradation_D-Xylose                      D-Xylose        
 9 D0308 0.00247    0.0438 Degradation  Sugar degradation_L-Rhamnose                    L-Rhamnose      
10 D0601 0.00370    0.0478 Degradation  Nitrogen compound degradation_Nitrate           Nitrate         
11 D0708 0.000165   0.0117 Degradation  Alcohol degradation_Phytol                      Phytol          

CI from acclimation to FMT1

element_gift %>%
  filter(time_point %in% c("FMT1","Acclimation") & type == "Treatment") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 2 × 6
  trait   p_value p_adjust Domain       Function                                Element           
  <chr>     <dbl>    <dbl> <chr>        <chr>                                   <chr>             
1 B0710 0.0000823   0.0119 Biosynthesis Vitamin biosynthesis_Phylloquinone (K1) Phylloquinone (K1)
2 D0102 0.000576    0.0418 Degradation  Lipid degradation_Fatty acid            Fatty acid        

10.1.3 CI vs WC

10.1.3.1 Alpha

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>%
  filter(type!="Control") %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point,  color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=8))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>%
  filter(type!="Control") 

Richness

Modelq0GLMMNB <- MASS::glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)

Call:
MASS::glm.nb(formula = richness ~ type * time_point, data = alpha_div_meta, 
    trace = TRUE, init.theta = 5.042866344, link = log)

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    4.4697     0.1527  29.279  < 2e-16 ***
typeTreatment                 -0.6557     0.2186  -2.999  0.00271 ** 
time_pointFMT1                -0.4209     0.2174  -1.936  0.05292 .  
typeTreatment:time_pointFMT1   0.5484     0.3146   1.743  0.08131 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(5.0429) family taken to be 1)

    Null deviance: 48.066  on 34  degrees of freedom
Residual deviance: 37.426  on 31  degrees of freedom
AIC: 336.8

Number of Fisher Scoring iterations: 1

              Theta:  5.04 
          Std. Err.:  1.34 

 2 x log-likelihood:  -326.796 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
 type        emmean    SE  df asymp.LCL asymp.UCL
 Hot_control   4.26 0.109 Inf      4.05      4.47
 Treatment     3.88 0.114 Inf      3.65      4.10

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE  df z.ratio p.value
 Hot_control - Treatment    0.381 0.157 Inf   2.425  0.0153

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 

Neutral

Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta) 
summary(Modelq1n)

Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.6988  -6.3363  -0.4774   7.6008  15.8281 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    44.601      3.324  13.418 1.88e-14 ***
typeTreatment                 -27.202      4.701  -5.786 2.26e-06 ***
time_pointFMT1                -18.862      4.701  -4.012 0.000353 ***
typeTreatment:time_pointFMT1   26.326      6.751   3.899 0.000483 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.972 on 31 degrees of freedom
Multiple R-squared:  0.5397,    Adjusted R-squared:  0.4951 
F-statistic: 12.12 on 3 and 31 DF,  p-value: 2.052e-05

Phylogenetic

Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta) 
summary(Model_phylo)

Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)

Residuals:
     Min       1Q   Median       3Q      Max 
-3.07987 -0.70068 -0.07723  0.54734  2.84255 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    6.5226     0.4274  15.260  5.8e-16 ***
typeTreatment                 -0.9826     0.6045  -1.626 0.114169    
time_pointFMT1                -2.2480     0.6045  -3.719 0.000793 ***
typeTreatment:time_pointFMT1   0.8981     0.8681   1.035 0.308891    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.282 on 31 degrees of freedom
Multiple R-squared:   0.39, Adjusted R-squared:  0.331 
F-statistic: 6.606 on 3 and 31 DF,  p-value: 0.001398

10.1.3.2 Beta diversity

Number of samples used

samples_to_keep_post4 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post4 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Control")
subset_meta_post4$type_time<-interaction(subset_meta_post4$type, subset_meta_post4$time_point)
length(samples_to_keep_post4)
[1] 35

Richness

richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
               colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     2 0.05438 0.027192 2.1499    999  0.126
Residuals 76 0.96123 0.012648                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
<<<<<<< HEAD
             Control Hot_control Treatment
Control                 0.718000     0.160
Hot_control 0.697512                 0.084
Treatment   0.155309    0.083221          
adonis2(phylo_post7 ~ type*time_point,
        data = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))),
        permutations = 999,
        strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))) %>% pull(individual),
          by="terms") %>%        
        broom::tidy() %>%
        tt()
======= Hot_control Treatment Hot_control 0.001 Treatment 0.00010134
adonis2(richness_post4 ~ type*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_zyvbojr9jfe6h7isy3qs, .table th.tinytable_css_zyvbojr9jfe6h7isy3qs { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_q3vlmkcbossuw9eaeyjx, .table th.tinytable_css_q3vlmkcbossuw9eaeyjx { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 2 0.2941494 0.07500331 3.592254 0.001
time_point 2 0.5528857 0.14097688 6.752032 0.001
type:time_point 4 0.2088312 0.05324856 1.275159 0.152
Residual 70 2.8659522 0.73077126 NA NA
Total 78 3.9218184 1.00000000 NA NA
<<<<<<< HEAD
subset_meta_post7_arrange <- column_to_rownames(subset_meta_post7, "Tube_code")
subset_meta_post7_arrange<-subset_meta_post7_arrange[labels(phylo_post7),]
pairwise <- pairwise.adonis(phylo_post7, subset_meta_post7_arrange$type, perm=999)
pairwise%>%
  tt()
=======
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(richness_post4),]

pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_xdkcihx41r32vi6u5eip, .table th.tinytable_css_xdkcihx41r32vi6u5eip { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_4l2ubp19yc5gisw70e31, .table th.tinytable_css_4l2ubp19yc5gisw70e31 { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs TreatmentTreatment.Acclimation vs Hot_control.Acclimation 1 1.3606630 5.087152 0.2412441 0.002 0.012 .
Treatment.Acclimation vs Treatment.FMT1 1 0.9551308 2.926054 0.1632291 0.001 0.006 *
Treatment.Acclimation vs Hot_control.FMT1 1 1.2263345 4.039487 0.2015764 0.001 0.006 *
Hot_control.Acclimation vs Treatment.FMT11 0.03324557 0.6023122 0.01190286 0.649 1.000
Control vs Hot_control 1 0.22899071 6.3999219 0.11149705 0.001 0.003 *
Treatment vs Hot_control 10.17684005 3.3769360 0.06210236 0.014 0.042 .0.3734921 1.268929 0.0779971 0.110 0.660
<<<<<<< HEAD
pairwise <- pairwise.adonis(phylo_post7, subset_meta_post7_arrange$time_point, perm=999)
pairwise%>%
  tt()
=======

Neutral

neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
               colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.087335 0.087335 15.924    999  0.003 **
Residuals 33 0.180987 0.005484                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.003
Treatment    0.00034565          
adonis2(neutral_post4 ~ type,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_xxa6ri3auop2jfg0ti55, .table th.tinytable_css_xxa6ri3auop2jfg0ti55 { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_vwsupixy2hc38u228396, .table th.tinytable_css_vwsupixy2hc38u228396 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 10.3794332 7.332180 0.12788943 0.001 0.003 *1.332543 0.1214599 4.562313 1
Acclimation vs FMT2 1 0.3351380 7.390585 0.12657152 0.001 0.003 *
FMT1 vs FMT2 1 0.1178756 3.274037 0.06032418 0.005 0.015 .
<<<<<<< HEAD

dbRDA

#Richness
cca_ord <- capscale(formula = richness_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post7, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post7$time_pointFMT1" = "FMT1",
                                 "subset_meta_post7$time_pointFMT2"="FMT2",
                                 "subset_meta_post7$typeHot_control"="Warm-control",
                                 "subset_meta_post7$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
                                 "subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
                                 "subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
                                 "subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")

beta_richness_nmds_post7 <-CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post7, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post7$time_pointFMT1" = "FMT1",
                                 "subset_meta_post7$time_pointFMT2"="FMT2",
                                 "subset_meta_post7$typeHot_control"="Warm-control",
                                 "subset_meta_post7$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
                                 "subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
                                 "subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
                                 "subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")

beta_neutral_nmds_post7 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()


#Phylogenetic
cca_ord <- capscale(formula = phylo_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post7, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post7$time_pointFMT1" = "FMT1",
                                 "subset_meta_post7$time_pointFMT2"="FMT2",
                                 "subset_meta_post7$typeHot_control"="Warm-control",
                                 "subset_meta_post7$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
                                 "subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
                                 "subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
                                 "subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")

beta_phylogenetic_nmds_post7 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()
ggarrange(beta_richness_nmds_post7, beta_neutral_nmds_post7, beta_phylogenetic_nmds_post7, ncol=3, nrow=1, common.legend = TRUE, legend="right")

10.2 What is the effect of FMT and temperature on the GM after 1 week compared to acclimation?

10.2.1 CI vs CC

10.2.1.1 Alpha diversity

label_map <- c(
  "Control" = "Cold-control", 
  "Treatment" = "Cold-intervention",
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity"
)
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(type %in% c("Control","Treatment") & time_point %in% c("FMT1","Acclimation")) %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
#  facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=10))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>%
  filter(type!="Hot_control") 

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(2.5327)  ( log )
Formula: richness ~ type * time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   329.2    338.4   -158.6    317.2       28 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4046 -0.7191  0.1331  0.6747  1.4330 

Random effects:
 Groups     Name        Variance  Std.Dev. 
 individual (Intercept) 2.062e-11 4.541e-06
Number of obs: 34, groups:  individual, 18

Fixed effects:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    3.9512     0.2275  17.368   <2e-16 ***
typeTreatment                 -0.1372     0.3132  -0.438    0.661    
time_pointFMT1                -0.3225     0.3140  -1.027    0.304    
typeTreatment:time_pointFMT1   0.4500     0.4435   1.015    0.310    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) typTrt t_FMT1
typeTretmnt -0.726              
tim_pntFMT1 -0.725  0.526       
typTr:_FMT1  0.513 -0.706 -0.708
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
 type      emmean    SE  df asymp.LCL asymp.UCL
 Control     3.79 0.157 Inf      3.48      4.10
 Treatment   3.88 0.157 Inf      3.57      4.18

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast            estimate    SE  df z.ratio p.value
 Control - Treatment  -0.0878 0.222 Inf  -0.396  0.6921

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE  df asymp.LCL asymp.UCL
 Acclimation   3.88 0.157 Inf      3.58      4.19
 FMT1          3.79 0.157 Inf      3.48      4.09

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE  df z.ratio p.value
 Acclimation - FMT1   0.0975 0.222 Inf   0.440  0.6603

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 

Neutral

Modelq1n <- lme(neutral ~ type*time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Modelq1n)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC   BIC    logLik
  244.9928 253.4 -116.4964

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    6.116057 8.521967

Fixed effects:  neutral ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                  21.027978  3.684717 16  5.706810  0.0000
typeTreatment                -3.628995  5.079636 16 -0.714420  0.4853
time_pointFMT1               -3.668645  4.182131 14 -0.877219  0.3952
typeTreatment:time_pointFMT1 10.917012  5.914427 14  1.845828  0.0862
 Correlation: 
                             (Intr) typTrt t_FMT1
typeTreatment                -0.725              
time_pointFMT1               -0.611  0.443       
typeTreatment:time_pointFMT1  0.432 -0.582 -0.707

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.79785802 -0.67295800  0.08387402  0.71606547  1.30935708 

Number of Observations: 34
Number of Groups: 18 
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
 type      emmean   SE df lower.CL upper.CL
 Control     19.2 2.92 14     12.9     25.5
 Treatment   21.0 2.92 14     14.8     27.3

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast            estimate   SE df t.ratio p.value
 Control - Treatment    -1.83 4.13 14  -0.443  0.6646

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation   19.2 2.54 16     13.8     24.6
 FMT1          21.0 2.54 14     15.6     26.5

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Acclimation - FMT1    -1.79 2.96 14  -0.605  0.5547

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

Phylogenetic

Model_phylo <- lme(phylogenetic ~ type*time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC    logLik
  125.282 133.6892 -56.64101

Random effects:
 Formula: ~1 | individual
         (Intercept) Residual
StdDev: 7.854299e-05 1.386163

Fixed effects:  phylogenetic ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   5.227341 0.4900827 16 10.666243  0.0000
typeTreatment                 0.312600 0.6735542 16  0.464105  0.6488
time_pointFMT1               -0.808097 0.6735542 14 -1.199751  0.2501
typeTreatment:time_pointFMT1 -0.541741 0.9525495 14 -0.568728  0.5786
 Correlation: 
                             (Intr) typTrt t_FMT1
typeTreatment                -0.728              
time_pointFMT1               -0.728  0.529       
typeTreatment:time_pointFMT1  0.514 -0.707 -0.707

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-2.630942795 -0.368275262  0.004976243  0.364227875  2.050663052 

Number of Observations: 34
Number of Groups: 18 
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
 type      emmean    SE df lower.CL upper.CL
 Control     4.82 0.337 14     4.10     5.55
 Treatment   4.87 0.337 14     4.14     5.59

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast            estimate    SE df t.ratio p.value
 Control - Treatment  -0.0417 0.476 14  -0.088  0.9314

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE df lower.CL upper.CL
 Acclimation   5.38 0.337 16     4.67     6.10
 FMT1          4.30 0.337 14     3.58     5.03

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE df t.ratio p.value
 Acclimation - FMT1     1.08 0.476 14   2.265  0.0399

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

10.2.1.2 Beta diversity

Number of samples used

samples_to_keep_post3 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post3 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control")
subset_meta_post3$time_point<-as.factor(subset_meta_post3$time_point)
subset_meta_post3$type<-as.factor(subset_meta_post3$type)
length(samples_to_keep_post3)
[1] 34

Richness

richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
               colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.002424 0.0024241 0.4717    999  0.482
Residuals 32 0.164445 0.0051389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Control Treatment
Control                0.47
Treatment 0.49714          
adonis2(richness_post3 ~ type*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
=======
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(neutral_post4),]
pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_xequglxj0vva2o6kdxtz, .table th.tinytable_css_xequglxj0vva2o6kdxtz { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_hxyn4amlq8uzaste8w29, .table th.tinytable_css_hxyn4amlq8uzaste8w29 { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 1 0.4269704 0.03612811 1.312996 0.001
time_point 1 1.1427754 0.09669595 3.514201 0.001
type:time_point 1 0.4928541 0.04170286 1.515598 0.050
Residual 30 9.7556341 0.82547308 NA NA
Total 33 11.8182341 1.00000000 NA NATreatment.FMT1 vs Hot_control.FMT1 1 0.4150076 1.637268 0.09840968 0.062 0.372
<<<<<<< HEAD
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(richness_post3),]

pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$type, perm=999)
pairwise%>%
  tt()
=======

Phylogenetic

phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
               colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07659 0.076588 5.0159    999  0.026 *
Residuals 33 0.50388 0.015269                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.033
Treatment      0.031971          
adonis2(phylo_post4 ~ type,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_l5js7rbeqaxhumvczvgh, .table th.tinytable_css_l5js7rbeqaxhumvczvgh { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_k4fj9c6v68b3z2lu952e, .table th.tinytable_css_k4fj9c6v68b3z2lu952e { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 10.4269704 1.199433 0.03612811 0.166 0.166 0.2070779 0.09373244 3.413088 1
Residual 33 2.0021670 0.90626756 NA NA
Total 34 2.2092449 1.00000000 NA NA
<<<<<<< HEAD
pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$time_point, perm=999)
pairwise%>%
  tt()
=======

NMDS

beta_richness_nmds_post4 <- richness_post4 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post4, by = join_by(sample == Tube_code)) %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_neutral_nmds_post4 <- neutral_post4 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post4, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_phylogenetic_nmds_post4 <- phylo_post4 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post4, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")
ggarrange(beta_richness_nmds_post4, beta_neutral_nmds_post4, beta_phylogenetic_nmds_post4, ncol=2, nrow=2, common.legend = TRUE, legend="right")

10.1.3.3 Functional differences

WC from acclimation to FMT1

element_gift %>%
  filter(time_point %in% c("FMT1","Acclimation") & type == "Hot_control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>

10.2 What is the effect of FMT and temperature on the GM after 2 weeks caompared to acclimation?

10.2.1 CI vs CC

10.2.1.1 Alpha

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(type %in% c("Control","Treatment") & time_point %in% c("FMT2","Acclimation")) %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
#  facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=10))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>%
  filter(type!="Hot_control") 

Richness

Modelq0GLMMNB <- MASS::glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)

Call:
MASS::glm.nb(formula = richness ~ type * time_point, data = alpha_div_meta, 
    trace = TRUE, init.theta = 4.087439892, link = log)

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    3.9512     0.1816  21.756   <2e-16 ***
typeTreatment                 -0.1372     0.2502  -0.548    0.584    
time_pointFMT2                 0.1206     0.2491   0.484    0.628    
typeTreatment:time_pointFMT2   0.4149     0.3469   1.196    0.232    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(4.0874) family taken to be 1)

    Null deviance: 43.336  on 34  degrees of freedom
Residual deviance: 37.889  on 31  degrees of freedom
AIC: 340.92

Number of Fisher Scoring iterations: 1

              Theta:  4.09 
          Std. Err.:  1.07 

 2 x log-likelihood:  -330.917 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
 type      emmean    SE  df asymp.LCL asymp.UCL
 Control     4.01 0.125 Inf      3.77      4.26
 Treatment   4.08 0.121 Inf      3.85      4.32

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast            estimate    SE  df z.ratio p.value
 Control - Treatment  -0.0702 0.173 Inf  -0.405  0.6855

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 

Neutral

Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta) 
summary(Modelq1n)

Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)

Residuals:
    Min      1Q  Median      3Q     Max 
-20.377  -6.596   2.212   6.374  21.492 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    21.373      3.672   5.820 2.05e-06 ***
typeTreatment                  -3.974      5.047  -0.787  0.43704    
time_pointFMT2                  2.345      5.047   0.465  0.64552    
typeTreatment:time_pointFMT2   20.138      7.032   2.864  0.00745 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.39 on 31 degrees of freedom
Multiple R-squared:  0.4388,    Adjusted R-squared:  0.3845 
F-statistic:  8.08 on 3 and 31 DF,  p-value: 0.0004037

Phylogenetic

Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta) 
summary(Model_phylo)

Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6469 -0.6788  0.0109  0.4658  2.9735 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   5.22734    0.52818   9.897  4.1e-11 ***
typeTreatment                 0.31260    0.72591   0.431    0.670    
time_pointFMT2               -0.45260    0.72591  -0.623    0.538    
typeTreatment:time_pointFMT2  0.02719    1.01138   0.027    0.979    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.494 on 31 degrees of freedom
Multiple R-squared:  0.03742,   Adjusted R-squared:  -0.05573 
F-statistic: 0.4017 on 3 and 31 DF,  p-value: 0.7527

10.2.1.2 Beta diversity

Number of samples used

samples_to_keep_post5 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post5 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control")
subset_meta_post5$type_time<-interaction(subset_meta_post5$type, subset_meta_post5$time_point)
length(samples_to_keep_post5)
[1] 35

Richness

richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
               colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.004294 0.0042942 0.4761    999  0.488
Residuals 33 0.297625 0.0090189                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Control Treatment
Control               0.486
Treatment 0.49501          
adonis2(richness_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_qnbqx8co8bmcd7x3kkzt, .table th.tinytable_css_qnbqx8co8bmcd7x3kkzt { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_fvxtb9e31nid36v4k606, .table th.tinytable_css_fvxtb9e31nid36v4k606 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 1 1.135334 3.400826 0.09606629 0.001 0.001 **

Neutral

neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
               colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.019686 0.0196859 2.8068    999  0.109
Residuals 32 0.224435 0.0070136                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Control Treatment
Control                 0.1
Treatment 0.10361          
adonis2(neutral_post3 ~ type,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual),
        by="terms") %>%        
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
type 1 0.4361022 0.04050485 1.350872 1
Residual 32 10.3305628 0.95949515 NA NA
Total 33 10.7666650 1.00000000 NA NA
<<<<<<< HEAD
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(neutral_post3),]
pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$type, perm=999)
pairwise%>%
  tt()
=======
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(richness_post5),]

pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_upmfaiopedcmahslh2sq, .table th.tinytable_css_upmfaiopedcmahslh2sq { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_grdma35xmar66z9luemd, .table th.tinytable_css_grdma35xmar66z9luemd { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs TreatmentControl.Acclimation vs Treatment.Acclimation 1 0.3657243 1.123239 0.06966584 0.250 1.000
Control.Acclimation vs Control.FMT21 0.4361022 1.350872 0.04050485 0.143 0.143
pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$time_point, perm=999)
pairwise%>%
  tt()
<<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1Control.Acclimation vs Treatment.FMT2 1 0.8072604 2.940901 0.16392161 0.003 0.018 .
Treatment.Acclimation vs Control.FMT2 1 0.5959230 1.984036 0.11032208 0.002 0.012 .
Treatment.Acclimation vs Treatment.FMT21 1.683105 5.929324 0.1563256 0.001 0.001 **
<<<<<<< HEAD

Phylogenetic

phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
               colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
=======

Neutral

neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
               colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
<<<<<<< HEAD
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.04697 0.046974 2.8076    999  0.112
Residuals 32 0.53540 0.016731                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Control Treatment
Control               0.123
Treatment 0.10357          
adonis2(phylo_post3 ~ type,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual),
        by="terms") %>%        
        broom::tidy() %>%
        tt()
======= Df Sum Sq Mean Sq F N.Perm Pr(>F) Groups 1 0.034255 0.034255 3.9748 999 0.055 . Residuals 33 0.284399 0.008618 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Pairwise comparisons: (Observed p-value below diagonal, permuted p-value above diagonal) Control Treatment Control 0.051 Treatment 0.054506
adonis2(neutral_post5 ~ type,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_n02xo7zp855yrblo2als, .table th.tinytable_css_n02xo7zp855yrblo2als { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_41s8p3wwr78e837y3ifn, .table th.tinytable_css_41s8p3wwr78e837y3ifn { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 1 0.03710337 0.01720754 0.5602822 10.642454 0.06322634 2.227293 1
Residual 32 2.11912470 0.98279246 NA NA
Total 33 2.15622808 1.00000000 NA NA
<<<<<<< HEAD
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(phylo_post3),]
pairwise <- pairwise.adonis(phylo_post3, subset_meta_post3_arrange$type, perm=999)
pairwise%>%
  tt()
=======
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(neutral_post5),]
pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_z3btk0dexo1og70pzmf6, .table th.tinytable_css_z3btk0dexo1og70pzmf6 { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_ao0irp265zjiku95dnnu, .table th.tinytable_css_ao0irp265zjiku95dnnu { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 10.03710337 0.5602822 0.01720754 0.678 0.6780.2759858 0.9928976 0.06208366 0.450 1.000
pairwise <- pairwise.adonis(phylo_post3, subset_meta_post3_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 1 0.3325691 5.835636 0.1542365 0.001 0.001 **

dbRDA

#Richness
cca_ord <- capscale(formula = richness_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post3, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post3$time_pointFMT1" = "FMT1",
                                 "subset_meta_post3$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")

beta_richness_nmds_post3 <-CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Treatment"),
                     labels=c("Cold-Cold", "Cold-Hot"),
                     values=c("#4477AA","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post3, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post3$time_pointFMT1" = "FMT1",
                                 "subset_meta_post3$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")

beta_neutral_nmds_post3 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Treatment"),
                     labels=c("Cold-Cold", "Cold-Hot"),
                     values=c("#4477AA","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()


#Phylogenetic
cca_ord <- capscale(formula = phylo_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post3, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post3$time_pointFMT1" = "FMT1",
                                 "subset_meta_post3$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")

beta_phylogenetic_nmds_post3 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Treatment"),
                     labels=c("Cold-Cold", "Cold-Hot"),
                     values=c("#4477AA","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()
ggarrange(beta_richness_nmds_post3, beta_neutral_nmds_post3, beta_phylogenetic_nmds_post3, ncol=3, nrow=1, common.legend = TRUE, legend="right")

10.2.2 Functional differences

CC from acclimation to FMT1

element_gift %>%
  filter(time_point %in% c("FMT1","Acclimation") & type == "Control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 11 × 6
   trait  p_value p_adjust Domain       Function                                        Element         
   <chr>    <dbl>    <dbl> <chr>        <chr>                                           <chr>           
 1 B0204 0.00370    0.0478 Biosynthesis Amino acid biosynthesis_Serine                  Serine          
 2 B0205 0.000165   0.0117 Biosynthesis Amino acid biosynthesis_Threonine               Threonine       
 3 B0703 0.00370    0.0478 Biosynthesis Vitamin biosynthesis_Niacin (B3)                Niacin (B3)     
 4 B0805 0.00156    0.0370 Biosynthesis Aromatic compound biosynthesis_Indole-3-acetate Indole-3-acetate
 5 D0102 0.000329   0.0156 Degradation  Lipid degradation_Fatty acid                    Fatty acid      
 6 D0212 0.000987   0.0350 Degradation  Polysaccharide degradation_Arabinan             Arabinan        
 7 D0304 0.00247    0.0438 Degradation  Sugar degradation_D-Arabinose                   D-Arabinose     
 8 D0306 0.00156    0.0370 Degradation  Sugar degradation_D-Xylose                      D-Xylose        
 9 D0308 0.00247    0.0438 Degradation  Sugar degradation_L-Rhamnose                    L-Rhamnose      
10 D0601 0.00370    0.0478 Degradation  Nitrogen compound degradation_Nitrate           Nitrate         
11 D0708 0.000165   0.0117 Degradation  Alcohol degradation_Phytol                      Phytol          

CI from acclimation to FMT1

element_gift %>%
  filter(time_point %in% c("FMT1","Acclimation") & type == "Treatment") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 2 × 6
  trait   p_value p_adjust Domain       Function                                Element           
  <chr>     <dbl>    <dbl> <chr>        <chr>                                   <chr>             
1 B0710 0.0000823   0.0119 Biosynthesis Vitamin biosynthesis_Phylloquinone (K1) Phylloquinone (K1)
2 D0102 0.000576    0.0418 Degradation  Lipid degradation_Fatty acid            Fatty acid        

10.2.3 CI vs WC

10.2.3.1 Alpha

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>%
  filter(type!="Control") %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point,  color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=8))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>%
  filter(type!="Control") 

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(5.0429)  ( log )
Formula: richness ~ type * time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   338.8    348.1   -163.4    326.8       29 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.85185 -0.65278  0.06667  0.54823  1.97668 

Random effects:
 Groups     Name        Variance  Std.Dev. 
 individual (Intercept) 1.315e-11 3.627e-06
Number of obs: 35, groups:  individual, 18

Fixed effects:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    4.4697     0.1527  29.279  < 2e-16 ***
typeTreatment                 -0.6557     0.2186  -2.999  0.00271 ** 
time_pointFMT1                -0.4208     0.2174  -1.936  0.05292 .  
typeTreatment:time_pointFMT1   0.5484     0.3146   1.743  0.08131 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) typTrt t_FMT1
typeTretmnt -0.698              
tim_pntFMT1 -0.702  0.490       
typTr:_FMT1  0.485 -0.695 -0.691
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
 type        emmean    SE  df asymp.LCL asymp.UCL
 Hot_control   4.26 0.109 Inf      4.05      4.47
 Treatment     3.88 0.114 Inf      3.65      4.10

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE  df z.ratio p.value
 Hot_control - Treatment    0.381 0.157 Inf   2.425  0.0153

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE  df asymp.LCL asymp.UCL
 Acclimation   4.14 0.109 Inf      3.93      4.36
 FMT1          4.00 0.113 Inf      3.77      4.22

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE  df z.ratio p.value
 Acclimation - FMT1    0.147 0.157 Inf   0.932  0.3512

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 

Neutral

Modelq1n <- lme(neutral ~ type*time_point, 
               random = ~ 1 | individual,
               data = alpha_div_meta,) 
summary(Modelq1n)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  250.0982 258.7021 -119.0491

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    5.122922 8.546652

Fixed effects:  neutral ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   44.60136  3.321472 16 13.428192   0e+00
typeTreatment                -27.20238  4.697271 16 -5.791103   0e+00
time_pointFMT1               -18.86154  4.028930 15 -4.681526   3e-04
typeTreatment:time_pointFMT1  26.15805  5.809237 15  4.502838   4e-04
 Correlation: 
                             (Intr) typTrt t_FMT1
typeTreatment                -0.707              
time_pointFMT1               -0.606  0.429       
typeTreatment:time_pointFMT1  0.421 -0.595 -0.694

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-1.903393212 -0.678106330  0.002306491  0.652929042  1.401732098 

Number of Observations: 35
Number of Groups: 18 
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
 type        emmean   SE df lower.CL upper.CL
 Hot_control   35.2 2.64 15     29.5     40.8
 Treatment     21.0 2.70 15     15.3     26.8

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                estimate   SE df t.ratio p.value
 Hot_control - Treatment     14.1 3.78 15   3.739  0.0020

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation   31.0 2.35 16     26.0     36.0
 FMT1          25.2 2.42 15     20.1     30.4

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate  SE df t.ratio p.value
 Acclimation - FMT1     5.78 2.9 15   1.991  0.0650

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

Phylogenetic

Model_phylo <- lme(phylogenetic ~ type*time_point, 
               random = ~ 1 | individual,
               data = alpha_div_meta) 
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  124.0635 132.6674 -56.03174

Random effects:
 Formula: ~1 | individual
         (Intercept) Residual
StdDev: 0.0002455305 1.282332

Fixed effects:  phylogenetic ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   6.522583 0.4274440 16 15.259503  0.0000
typeTreatment                -0.982643 0.6044971 16 -1.625554  0.1236
time_pointFMT1               -2.247960 0.6044971 15 -3.718727  0.0021
typeTreatment:time_pointFMT1  0.898121 0.8681429 15  1.034531  0.3173
 Correlation: 
                             (Intr) typTrt t_FMT1
typeTreatment                -0.707              
time_pointFMT1               -0.707  0.500       
typeTreatment:time_pointFMT1  0.492 -0.696 -0.696

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.40177643 -0.54641173 -0.06022755  0.42683546  2.21670618 

Number of Observations: 35
Number of Groups: 18 
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
 type        emmean    SE df lower.CL upper.CL
 Hot_control   5.40 0.302 15     4.75     6.04
 Treatment     4.87 0.312 15     4.20     5.53

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE df t.ratio p.value
 Hot_control - Treatment    0.534 0.434 15   1.229  0.2379

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE df lower.CL upper.CL
 Acclimation   6.03 0.302 16     5.39     6.67
 FMT1          4.23 0.312 15     3.57     4.90

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE df t.ratio p.value
 Acclimation - FMT1      1.8 0.434 15   4.144  0.0009

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

10.2.3.2 Beta diversity

Number of samples used

samples_to_keep_post4 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post4 <- sample_metadata %>%
  filter(time_point=="FMT1"|time_point=="Acclimation") %>% 
  filter(type!="Control")
subset_meta_post4$time_point<-as.factor(subset_meta_post4$time_point)
subset_meta_post4$type<-as.factor(subset_meta_post4$type)
length(samples_to_keep_post4)
[1] 35

Richness

richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
               colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.076276 0.076276 19.518    999  0.001 ***
Residuals 33 0.128967 0.003908                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.001
Treatment    0.00010134          
adonis2(richness_post4 ~ type*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
type 1 1.1360754 0.09995918 4.050611 0.001
time_point 1 0.9251516 0.08140076 3.298575 0.001
type:time_point 1 0.6095926 0.05363586 2.173467 0.002
Residual 31 8.6945735 0.76500420 NA NA
Total 34 11.3653932 1.00000000 NA NA
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(richness_post4),]

pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Hot_control 1 1.136075 3.665004 0.09995918 0.001 0.001 **
pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 1 0.9366646 2.963922 0.08241375 0.001 0.001 **

Neutral

neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
               colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.087335 0.087335 15.924    999  0.002 **
Residuals 33 0.180987 0.005484                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.001
Treatment    0.00034565          
adonis2(neutral_post4 ~ type,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
type 1 1.332543 0.1214599 4.562313 1
Residual 33 9.638511 0.8785401 NA NA
Total 34 10.971054 1.0000000 NA NA
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(neutral_post4),]
pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Hot_control 1 1.332543 4.562313 0.1214599 0.001 0.001 **
pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 1 1.291025 4.401211 0.1176756 0.001 0.001 **

Phylogenetic

phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
               colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07659 0.076588 5.0159    999  0.035 *
Residuals 33 0.50388 0.015269                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.036
Treatment      0.031971          
adonis2(phylo_post4 ~ type,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
type 1 0.2070779 0.09373244 3.413088 1
Residual 33 2.0021670 0.90626756 NA NA
Total 34 2.2092449 1.00000000 NA NA
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(phylo_post4),]
pairwise <- pairwise.adonis(phylo_post4, subset_meta_post4_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Hot_control 1 0.2070779 3.413088 0.09373244 0.012 0.012 .
pairwise <- pairwise.adonis(phylo_post4, subset_meta_post4_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT1 1 0.2895801 4.978027 0.1310765 0.002 0.002 *

RDA

#Richness
cca_ord <- capscale(formula = richness_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post4, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post4$time_pointFMT1" = "FMT1",
                                 "subset_meta_post4$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")

beta_richness_nmds_post4 <-CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post4, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post4$time_pointFMT1" = "FMT1",
                                 "subset_meta_post4$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")

beta_neutral_nmds_post4 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()


#Phylogenetic
cca_ord <- capscale(formula = phylo_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post4, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post4$time_pointFMT1" = "FMT1",
                                 "subset_meta_post4$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")

beta_phylogenetic_nmds_post4 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()
ggarrange(beta_richness_nmds_post4, beta_neutral_nmds_post4, beta_phylogenetic_nmds_post4, ncol=3, nrow=1, common.legend = TRUE, legend="right")

10.2.3.3 Functional differences

WC from acclimation to FMT1

element_gift %>%
  filter(time_point %in% c("FMT1","Acclimation") & type == "Hot_control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>

10.3 What is the effect of FMT and temperature on the GM after 2 weeks caompared to acclimation?

10.3.1 CI vs CC

10.3.1.1 Alpha

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(type %in% c("Control","Treatment") & time_point %in% c("FMT2","Acclimation")) %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
#  facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control", "Treatment"),
          values=c("#4477AA50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=10))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>%
  filter(type!="Hot_control") 

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(4.0875)  ( log )
Formula: richness ~ type * time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   342.9    352.2   -165.5    330.9       29 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7595 -0.2361  0.2110  0.4983  1.1980 

Random effects:
 Groups     Name        Variance  Std.Dev. 
 individual (Intercept) 6.519e-12 2.553e-06
Number of obs: 35, groups:  individual, 18

Fixed effects:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    3.9512     0.1816  21.756   <2e-16 ***
typeTreatment                 -0.1372     0.2502  -0.548    0.583    
time_pointFMT2                 0.1206     0.2491   0.484    0.628    
typeTreatment:time_pointFMT2   0.4149     0.3469   1.196    0.232    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) typTrt t_FMT2
typeTretmnt -0.726              
tim_pntFMT2 -0.729  0.529       
typTr:_FMT2  0.524 -0.721 -0.718
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
 type      emmean    SE  df asymp.LCL asymp.UCL
 Control     4.01 0.125 Inf      3.77      4.26
 Treatment   4.08 0.121 Inf      3.85      4.32

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast            estimate    SE  df z.ratio p.value
 Control - Treatment  -0.0702 0.173 Inf  -0.405  0.6855

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE  df asymp.LCL asymp.UCL
 Acclimation   3.88 0.125 Inf      3.64      4.13
 FMT2          4.21 0.120 Inf      3.98      4.45

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE  df z.ratio p.value
 Acclimation - FMT2   -0.328 0.173 Inf  -1.892  0.0585

Results are averaged over the levels of: type 
Results are given on the log (not the response) scale. 

Neutral

Modelq1n <- lme(neutral ~ type*time_point, 
               random = ~ 1 | individual,
               data = alpha_div_meta)  
summary(Modelq1n)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  253.6069 262.2108 -120.8034

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.268167 9.858563

Fixed effects:  neutral ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                  21.455542  3.670059 16  5.846103  0.0000
typeTreatment                -4.056558  5.045308 16 -0.804026  0.4332
time_pointFMT2                2.262002  4.804331 15  0.470826  0.6445
typeTreatment:time_pointFMT2 20.220508  6.684284 15  3.025082  0.0085
 Correlation: 
                             (Intr) typTrt t_FMT2
typeTreatment                -0.727              
time_pointFMT2               -0.697  0.507       
typeTreatment:time_pointFMT2  0.501 -0.684 -0.719

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.7676771 -0.6504757  0.2922823  0.6462932  2.0120364 

Number of Observations: 35
Number of Groups: 18 
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
 type      emmean   SE df lower.CL upper.CL
 Control     22.6 2.64 15     17.0     28.2
 Treatment   28.6 2.57 15     23.2     34.1

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast            estimate   SE df t.ratio p.value
 Control - Treatment    -6.05 3.68 15  -1.645  0.1208

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation   19.4 2.52 16     14.1     24.8
 FMT2          31.8 2.45 15     26.6     37.0

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Acclimation - FMT2    -12.4 3.34 15  -3.702  0.0021

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

Phylogenetic

Model_phylo <- lme(phylogenetic ~ type*time_point, 
               random = ~ 1 | individual,
               data = alpha_div_meta)  
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC    logLik
  133.508 142.1119 -60.75399

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.2974732 1.463985

Fixed effects:  phylogenetic ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   5.233354 0.5281280 16  9.909252  0.0000
typeTreatment                 0.306587 0.7258724 16  0.422370  0.6784
time_pointFMT2               -0.458617 0.7121982 15 -0.643946  0.5293
typeTreatment:time_pointFMT2  0.033201 0.9917181 15  0.033479  0.9737
 Correlation: 
                             (Intr) typTrt t_FMT2
typeTreatment                -0.728              
time_pointFMT2               -0.715  0.521       
typeTreatment:time_pointFMT2  0.514 -0.705 -0.718

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-2.403145672 -0.437916442 -0.005178256  0.313089271  1.956678742 

Number of Observations: 35
Number of Groups: 18 
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
 type      emmean    SE df lower.CL upper.CL
 Control     5.00 0.370 15     4.22     5.79
 Treatment   5.33 0.359 15     4.56     6.09

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast            estimate    SE df t.ratio p.value
 Control - Treatment   -0.323 0.515 15  -0.627  0.5400

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE df lower.CL upper.CL
 Acclimation   5.39 0.363 16     4.62     6.16
 FMT2          4.94 0.352 15     4.19     5.70

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE df t.ratio p.value
 Acclimation - FMT2    0.442 0.496 15   0.891  0.3868

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

10.3.1.2 Beta diversity

Number of samples used

samples_to_keep_post5 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post5 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Hot_control")
subset_meta_post5$time_point<-as.factor(subset_meta_post5$time_point)
subset_meta_post5$type<-as.factor(subset_meta_post5$type)
length(samples_to_keep_post5)
[1] 35

Richness

richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
               colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.004294 0.0042942 0.4761    999   0.52
Residuals 33 0.297625 0.0090189                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Control Treatment
Control               0.524
Treatment 0.49501          
adonis2(richness_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
type 1 0.5174241 0.04775345 1.797587 0.001
time_point 1 0.9068994 0.08369841 3.150666 0.001
type:time_point 1 0.4878440 0.04502349 1.694822 0.028
Residual 31 8.9231562 0.82352465 NA NA
Total 34 10.8353238 1.00000000 NA NA
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(richness_post5),]

pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.5174241 1.65489 0.04775345 0.022 0.022 .
pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT2 1 0.9000618 2.989558 0.08306737 0.001 0.001 **

Neutral

neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
               colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.034255 0.034255 3.9748    999  0.063 .
Residuals 33 0.284399 0.008618                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Control Treatment
Control                0.058
Treatment 0.054506          
adonis2(neutral_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
type 1 0.6424540 0.06322634 2.589936 0.001
time_point 1 1.1894033 0.11705370 4.794863 0.001
type:time_point 1 0.6395258 0.06293816 2.578132 0.010
Residual 31 7.6897933 0.75678179 NA NA
Total 34 10.1611764 1.00000000 NA NA
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(neutral_post5),]
pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.642454 2.227293 0.06322634 0.007 0.007 *
pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT2 1 1.172649 4.305202 0.1154049 0.001 0.001 **

Phylogenetic

phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
               colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02670 0.026703 1.6775    999  0.185
Residuals 33 0.52532 0.015919                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Control Treatment
Control               0.205
Treatment 0.20425          
adonis2(phylo_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
type 1 0.03644585 0.01867131 0.7636208 0.003
time_point 1 0.35441680 0.18156873 7.4258113 0.002
type:time_point 1 0.08154954 0.04177806 1.7086422 0.204
Residual 31 1.47955831 0.75798190 NA NA
Total 34 1.95197051 1.00000000 NA NA
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(phylo_post5),]
pairwise <- pairwise.adonis(phylo_post5, subset_meta_post5_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.03644585 0.6278766 0.01867131 0.603 0.603
pairwise <- pairwise.adonis(phylo_post5, subset_meta_post5_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT2 1 0.3562316 7.366895 0.1824984 0.001 0.001 **

dbRDA

#Richness
cca_ord <- capscale(formula = richness_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post5, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post5$time_pointFMT2" = "FMT2",
                                 "subset_meta_post5$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")

beta_richness_nmds_post5 <-CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Treatment"),
                     labels=c("Cold-Cold", "Cold-Hot"),
                     values=c("#4477AA","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post5, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post5$time_pointFMT2" = "FMT2",
                                 "subset_meta_post5$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")

beta_neutral_nmds_post5 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Treatment"),
                     labels=c("Cold-Cold", "Cold-Hot"),
                     values=c("#4477AA","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()


#Phylogenetic
cca_ord <- capscale(formula = phylo_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post5, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post5$time_pointFMT2" = "FMT2",
                                 "subset_meta_post5$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")

beta_phylogenetic_nmds_post5 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Treatment"),
                     labels=c("Cold-Cold", "Cold-Hot"),
                     values=c("#4477AA","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()
ggarrange(beta_richness_nmds_post5, beta_neutral_nmds_post5, beta_phylogenetic_nmds_post5, ncol=3, nrow=1, common.legend = TRUE, legend="right")

10.3.1.3 Functional differences

CC from acclimation to FMT2

element_gift %>%
  filter(time_point %in% c("FMT2","Acclimation") & type == "Control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>

CI from acclimation to FMT2

element_gift %>%
  filter(time_point %in% c("FMT2","Acclimation") & type == "Treatment") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 17 × 6
   trait  p_value p_adjust Domain       Function                               Element           
   <chr>    <dbl>    <dbl> <chr>        <chr>                                  <chr>             
 1 B0103 0.000787   0.0230 Biosynthesis Nucleic acid biosynthesis_UDP/UTP      UDP/UTP           
 2 B0214 0.000494   0.0230 Biosynthesis Amino acid biosynthesis_Glutamate      Glutamate         
 3 D0102 0.000782   0.0230 Degradation  Lipid degradation_Fatty acid           Fatty acid        
 4 D0205 0.00399    0.0448 Degradation  Polysaccharide degradation_Pectin      Pectin            
 5 D0210 0.00123    0.0234 Degradation  Polysaccharide degradation_Beta-mannan Beta-mannan       
 6 D0308 0.00276    0.0366 Degradation  Sugar degradation_L-Rhamnose           L-Rhamnose        
 7 D0512 0.000494   0.0230 Degradation  Amino acid degradation_Histidine       Histidine         
 8 D0513 0.00564    0.0484 Degradation  Amino acid degradation_Tryptophan      Tryptophan        
 9 D0601 0.000494   0.0230 Degradation  Nitrogen compound degradation_Nitrate  Nitrate           
10 D0613 0.00564    0.0484 Degradation  Nitrogen compound degradation_Taurine  Taurine           
11 D0706 0.00564    0.0484 Degradation  Alcohol degradation_Ethylene glycol    Ethylene glycol   
12 D0708 0.00185    0.0270 Degradation  Alcohol degradation_Phytol             Phytol            
13 D0801 0.00145    0.0234 Degradation  Xenobiotic degradation_Toluene         Toluene           
14 D0802 0.00145    0.0234 Degradation  Xenobiotic degradation_Xylene          Xylene            
15 D0901 0.00564    0.0484 Degradation  Antibiotic degradation_Penicillin      Penicillin        
16 D0902 0.00108    0.0234 Degradation  Antibiotic degradation_Carbapenem      Carbapenem        
17 S0105 0.00399    0.0448 Structure    Cellular structure_Lipopolysaccharide  Lipopolysaccharide

10.3.2 CI vs WC

10.3.2.1 Alpha

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>%
  filter(type!="Control") %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point,  color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=8))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>%
  filter(type!="Control") 

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
emmeans(Modelq0GLMMNB, pairwise ~ type)
emmeans(Modelq0GLMMNB, pairwise ~ time_point)

Neutral

Modelq1n <- lme(neutral ~ type*time_point,
               random = ~ 1 | individual,
               data = alpha_div_meta)   
summary(Modelq1n)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  255.5098 264.3042 -121.7549

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    1.317597  9.38306

Fixed effects:  neutral ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   44.60136  3.158373 16 14.121626  0.0000
typeTreatment                -27.20238  4.466614 16 -6.090157  0.0000
time_pointFMT2               -12.40101  4.423217 16 -2.803618  0.0127
typeTreatment:time_pointFMT2  34.88352  6.255373 16  5.576569  0.0000
 Correlation: 
                             (Intr) typTrt t_FMT2
typeTreatment                -0.707              
time_pointFMT2               -0.700  0.495       
typeTreatment:time_pointFMT2  0.495 -0.700 -0.707

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.1054551 -0.5575213  0.2345834  0.5912103  2.0213400 

Number of Observations: 36
Number of Groups: 18 
emmeans(Modelq1n, pairwise ~ type)
$emmeans
 type        emmean   SE df lower.CL upper.CL
 Hot_control   38.4 2.25 16     33.6     43.2
 Treatment     28.6 2.25 16     23.9     33.4

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                estimate   SE df t.ratio p.value
 Hot_control - Treatment     9.76 3.19 16   3.061  0.0075

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation     31 2.23 16     26.3     35.7
 FMT2            36 2.23 16     31.3     40.8

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Acclimation - FMT2    -5.04 3.13 16  -1.612  0.1266

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

Phylogenetic

Model_phylo <- lme(phylogenetic ~ type*time_point, 
               random = ~ 1 | individual,
               data = alpha_div_meta)   
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  132.8692 141.6636 -60.43459

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.7604222 1.204042

Fixed effects:  phylogenetic ~ type * time_point 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   6.522583 0.4746882 16 13.740774  0.0000
typeTreatment                -0.982643 0.6713105 16 -1.463768  0.1626
time_pointFMT2               -1.066614 0.5675911 16 -1.879195  0.0786
typeTreatment:time_pointFMT2  0.641199 0.8026950 16  0.798807  0.4361
 Correlation: 
                             (Intr) typTrt t_FMT2
typeTreatment                -0.707              
time_pointFMT2               -0.598  0.423       
typeTreatment:time_pointFMT2  0.423 -0.598 -0.707

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.90560105 -0.54341938 -0.06251794  0.39498845  1.94321495 

Number of Observations: 36
Number of Groups: 18 
emmeans(Model_phylo, pairwise ~ type)
$emmeans
 type        emmean    SE df lower.CL upper.CL
 Hot_control   5.99 0.381 16     5.18     6.80
 Treatment     5.33 0.381 16     4.52     6.13

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE df t.ratio p.value
 Hot_control - Treatment    0.662 0.538 16   1.230  0.2364

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE df lower.CL upper.CL
 Acclimation   6.03 0.336 16     5.32     6.74
 FMT2          5.29 0.336 16     4.57     6.00

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE df t.ratio p.value
 Acclimation - FMT2    0.746 0.401 16   1.859  0.0816

Results are averaged over the levels of: type 
Degrees-of-freedom method: containment 

10.3.2.2 Beta diversity

Number of samples used

samples_to_keep_post6 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post6 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Control")
subset_meta_post6$time_point<-as.factor(subset_meta_post6$time_point)
subset_meta_post6$type<-as.factor(subset_meta_post6$type)
length(samples_to_keep_post6)
[1] 36

Richness

richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
               colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.054351 0.054351 9.1422    999  0.004 **
Residuals 34 0.202132 0.005945                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.006
Treatment     0.0047271          
adonis2(richness_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
<<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 11.2570893 0.11659822 4.856166 0.0010.4748999 1.9609749 0.11561688 0.002 0.012 .
time_point 1 0.6574224 0.06097760 2.539639 0.002
type:time_point 10.5831993 0.05409322 2.252913 0.004
Residual 32 8.2836650 0.76833097 NA NA
Total 35 10.7813760 1.00000000 NA NA0.6311089 2.4041625 0.13063146 0.002 0.012 .
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(richness_post6),]

pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Hot_control 1 1.257089 4.487584 0.1165982 0.001 0.001 **
pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT2 1 0.6574224 2.207869 0.0609776 0.003 0.003 *
<<<<<<< HEAD

Neutral

neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
               colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
=======

Phylogenetic

phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
               colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
<<<<<<< HEAD
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.078328 0.078328 12.006    999  0.001 ***
Residuals 34 0.221814 0.006524                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.005
Treatment     0.0014541          
adonis2(neutral_post6 ~ type,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
======= Df Sum Sq Mean Sq F N.Perm Pr(>F) Groups 1 0.02670 0.026703 1.6775 999 0.211 Residuals 33 0.52532 0.015919 Pairwise comparisons: (Observed p-value below diagonal, permuted p-value above diagonal) Control Treatment Control 0.233 Treatment 0.20425
adonis2(phylo_post5 ~ type,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_9qk78kjt0i9lfnp4ejpp, .table th.tinytable_css_9qk78kjt0i9lfnp4ejpp { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_25z8w88t8nonj4rkmnqn, .table th.tinytable_css_25z8w88t8nonj4rkmnqn { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 1 1.325563 0.1261942 4.910249 10.03644585 0.01867131 0.6278766 1
Residual 34 9.178585 0.8738058 NA NA
Total 35 10.504148 1.0000000 NA NA
<<<<<<< HEAD
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(neutral_post6),]
pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$type, perm=999)
pairwise%>%
  tt()
=======

NMDS

beta_richness_nmds_post5 <- richness_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code)) %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Treatment"),
                       labels=c("Cold-Cold", "Cold-Hot"),
                       values=c("#4477AA","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_neutral_nmds_post5 <- neutral_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Treatment"),
                       labels=c("Cold-Cold", "Cold-Hot"),
                       values=c("#4477AA","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_phylogenetic_nmds_post5 <- phylo_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Treatment"),
                       labels=c("Cold-Cold", "Cold-Hot"),
                       values=c("#4477AA","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")
ggarrange(beta_richness_nmds_post5, beta_neutral_nmds_post5, beta_phylogenetic_nmds_post5, ncol=2, nrow=2, common.legend = TRUE, legend="right")

10.2.1.3 Functional differences

CC from acclimation to FMT2

element_gift %>%
  filter(time_point %in% c("FMT2","Acclimation") & type == "Control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>

CI from acclimation to FMT2

element_gift %>%
  filter(time_point %in% c("FMT2","Acclimation") & type == "Treatment") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 17 × 6
   trait  p_value p_adjust Domain       Function                               Element           
   <chr>    <dbl>    <dbl> <chr>        <chr>                                  <chr>             
 1 B0103 0.000787   0.0230 Biosynthesis Nucleic acid biosynthesis_UDP/UTP      UDP/UTP           
 2 B0214 0.000494   0.0230 Biosynthesis Amino acid biosynthesis_Glutamate      Glutamate         
 3 D0102 0.000782   0.0230 Degradation  Lipid degradation_Fatty acid           Fatty acid        
 4 D0205 0.00399    0.0448 Degradation  Polysaccharide degradation_Pectin      Pectin            
 5 D0210 0.00123    0.0234 Degradation  Polysaccharide degradation_Beta-mannan Beta-mannan       
 6 D0308 0.00276    0.0366 Degradation  Sugar degradation_L-Rhamnose           L-Rhamnose        
 7 D0512 0.000494   0.0230 Degradation  Amino acid degradation_Histidine       Histidine         
 8 D0513 0.00564    0.0484 Degradation  Amino acid degradation_Tryptophan      Tryptophan        
 9 D0601 0.000494   0.0230 Degradation  Nitrogen compound degradation_Nitrate  Nitrate           
10 D0613 0.00564    0.0484 Degradation  Nitrogen compound degradation_Taurine  Taurine           
11 D0706 0.00564    0.0484 Degradation  Alcohol degradation_Ethylene glycol    Ethylene glycol   
12 D0708 0.00185    0.0270 Degradation  Alcohol degradation_Phytol             Phytol            
13 D0801 0.00145    0.0234 Degradation  Xenobiotic degradation_Toluene         Toluene           
14 D0802 0.00145    0.0234 Degradation  Xenobiotic degradation_Xylene          Xylene            
15 D0901 0.00564    0.0484 Degradation  Antibiotic degradation_Penicillin      Penicillin        
16 D0902 0.00108    0.0234 Degradation  Antibiotic degradation_Carbapenem      Carbapenem        
17 S0105 0.00399    0.0448 Structure    Cellular structure_Lipopolysaccharide  Lipopolysaccharide

10.2.2 CI vs WC

10.2.2.1 Alpha

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>%
  filter(type!="Control") %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
  ggplot(aes(y = value, x = time_point,  color = type, fill = type)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
  facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  scale_color_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Hot_control", "Treatment"),
          values=c("#d57d2c50","#76b18350")) +
  theme_minimal() +
  theme(axis.text.x=element_text(size=8))+
  labs(x = "Time Point", y = "value")

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>%
  filter(type!="Control") 

Richness

Modelq0GLMMNB <- glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)
emmeans(Modelq0GLMMNB, pairwise ~ type)

Neutral

Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta) 
summary(Modelq1n)

Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)

Residuals:
    Min      1Q  Median      3Q     Max 
-20.377  -5.358   2.158   5.663  19.246 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    44.601      3.158  14.122 2.68e-15 ***
typeTreatment                 -27.202      4.467  -6.090 8.36e-07 ***
time_pointFMT2                -12.401      4.467  -2.776  0.00911 ** 
typeTreatment:time_pointFMT2   34.884      6.317   5.522 4.34e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.475 on 32 degrees of freedom
Multiple R-squared:  0.571, Adjusted R-squared:  0.5308 
F-statistic:  14.2 on 3 and 32 DF,  p-value: 4.69e-06

Phylogenetic

Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta) 
summary(Model_phylo)

Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0799 -0.7676 -0.0971  0.5549  2.9735 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    6.5226     0.4747  13.741  5.7e-15 ***
typeTreatment                 -0.9826     0.6713  -1.464    0.153    
time_pointFMT2                -1.0666     0.6713  -1.589    0.122    
typeTreatment:time_pointFMT2   0.6412     0.9494   0.675    0.504    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.424 on 32 degrees of freedom
Multiple R-squared:  0.1321,    Adjusted R-squared:  0.05075 
F-statistic: 1.624 on 3 and 32 DF,  p-value: 0.2033

10.2.2.2 Beta diversity

Number of samples used

samples_to_keep_post6 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Control") %>%
  select(Tube_code) %>% 
  pull()
subset_meta_post6 <- sample_metadata %>%
  filter(time_point=="FMT2"|time_point=="Acclimation") %>% 
  filter(type!="Control")
subset_meta_post6$type_time<-interaction(subset_meta_post6$type, subset_meta_post6$time_point)
length(samples_to_keep_post6)
[1] 36

Richness

richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
               colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.054351 0.054351 9.1422    999  0.005 **
Residuals 34 0.202132 0.005945                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.005
Treatment     0.0047271          
adonis2(richness_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_f0zf6n7lz8158paups9l, .table th.tinytable_css_f0zf6n7lz8158paups9l { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_008xqo3uvzkbpi6zctrt, .table th.tinytable_css_008xqo3uvzkbpi6zctrt { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Hot_control 1 1.325563 4.910249 0.1261942 0.001 0.001 **
<<<<<<< HEAD
pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$time_point, perm=999)
pairwise%>%
  tt()
=======
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(richness_post6),]

pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_ruz9z9h43gtq848ow32t, .table th.tinytable_css_ruz9z9h43gtq848ow32t { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; } .table td.tinytable_css_ob11bhwyfi3ccm51mjew, .table th.tinytable_css_ob11bhwyfi3ccm51mjew { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT2Treatment.Acclimation vs Hot_control.Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.006 *
Treatment.Acclimation vs Treatment.FMT2 1 0.9130048 3.195028 0.16645080 0.001 0.006 *
Treatment.Acclimation vs Hot_control.FMT2 1 1.2747787 4.275366 0.21086503 0.001 0.006 *
Hot_control.Acclimation vs Treatment.FMT2 1 0.6397330 2.913695 0.15405213 0.001 0.006 *
Hot_control.Acclimation vs Hot_control.FMT2 1 0.3276169 1.412318 0.08111028 0.042 0.252
Treatment.FMT2 vs Hot_control.FMT21 0.9072453 3.214197 0.0863702 0.002 0.002 *
<<<<<<< HEAD

Phylogenetic

phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
               colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
=======

Neutral

neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
               colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
<<<<<<< HEAD
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.04938 0.049381 3.8729    999  0.056 .
Residuals 34 0.43351 0.012750                       
=======
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.078328 0.078328 12.006    999  0.002 **
Residuals 34 0.221814 0.006524                        
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
<<<<<<< HEAD
Hot_control                  0.04
Treatment      0.057271          
adonis2(phylo_post6 ~ type,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
======= Hot_control 0.004 Treatment 0.0014541
adonis2(neutral_post6 ~ type,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_ov8p6edmajbwrv2eq9zh, .table th.tinytable_css_ov8p6edmajbwrv2eq9zh { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_22rqyvvcvazqw37ciq4q, .table th.tinytable_css_22rqyvvcvazqw37ciq4q { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
term df SumOfSqs R2 statistic p.value
type 1 0.1600877 0.08395907 3.116245 11.325563 0.1261942 4.910249 1
Residual 34 1.7466474 0.91604093 NA NA
Total 35 1.9067351 1.00000000 NA NA
<<<<<<< HEAD
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(phylo_post6),]
pairwise <- pairwise.adonis(phylo_post6, subset_meta_post6_arrange$type, perm=999)
pairwise%>%
  tt()
=======
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(neutral_post6),]
pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_vj532nvb8indd1jitiw0, .table th.tinytable_css_vj532nvb8indd1jitiw0 { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_9qbyiug1baoth4sa8w1r, .table th.tinytable_css_9qbyiug1baoth4sa8w1r { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Hot_control 10.1600877 3.116245 0.08395907 0.016 0.016 .
1.6125755 5.982598 0.27215155 0.001 0.006 *
Hot_control.Acclimation vs Treatment.FMT2 1 0.6202327 3.151987 0.16457754 0.001 0.006 *
Hot_control.Acclimation vs Hot_control.FMT2 1 0.3634438 1.708339 0.09647087 0.026 0.156
Treatment.FMT2 vs Hot_control.FMT2 1 0.5010202 2.206532 0.12119453 0.001 0.006 *
<<<<<<< HEAD
pairwise <- pairwise.adonis(phylo_post6, subset_meta_post6_arrange$time_point, perm=999)
pairwise%>%
  tt()
=======

Phylogenetic

phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
               colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.04938 0.049381 3.8729    999  0.048 *
Residuals 34 0.43351 0.012750                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Hot_control Treatment
Hot_control                 0.034
Treatment      0.057271          
adonis2(phylo_post6 ~ type,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
======= .table td.tinytable_css_nfn92de6sdn79asszwrh, .table th.tinytable_css_nfn92de6sdn79asszwrh { border-bottom: solid #d3d8dc 0.1em; } .table td.tinytable_css_7xwp1sx2dsfixwp5fgpa, .table th.tinytable_css_7xwp1sx2dsfixwp5fgpa { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2 <<<<<<< HEAD ======= >>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs FMT2 10.2882491 6.055333 0.1511742 0.001 0.001 **0.1600877 0.08395907 3.116245 1
Residual 34 1.7466474 0.91604093 NA NA
Total 35 1.9067351 1.00000000 NA NA
<<<<<<< HEAD

dbRDA

#Richness
cca_ord <- capscale(formula = richness_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post6, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post6$time_pointFMT2" = "FMT2",
                                 "subset_meta_post6$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")

beta_richness_nmds_post6 <-CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post6, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post6$time_pointFMT2" = "FMT2",
                                 "subset_meta_post6$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")

beta_neutral_nmds_post6 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()


#Phylogenetic
cca_ord <- capscale(formula = phylo_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta_post6, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta_post6$time_pointFMT2" = "FMT2",
                                 "subset_meta_post6$typeTreatment" = "Cold-intervention",
                                 "subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")

beta_phylogenetic_nmds_post6 <- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
  scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()
ggarrange(beta_richness_nmds_post6, beta_neutral_nmds_post6, beta_phylogenetic_nmds_post6, ncol=3, nrow=1, common.legend = TRUE, legend="right")

=======

NMDS

beta_richness_nmds_post6 <- richness_post6 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post6, by = join_by(sample == Tube_code)) %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_neutral_nmds_post6 <- neutral_post6 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post6, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")

beta_phylogenetic_nmds_post6 <- phylo_post6 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post6, by = join_by(sample == Tube_code))%>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Hot_control", "Treatment"),
                       labels=c("Warm-control", "Cold-intervention"),
                       values=c("#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")
ggarrange(beta_richness_nmds_post6, beta_neutral_nmds_post6, beta_phylogenetic_nmds_post6, ncol=2, nrow=2, common.legend = TRUE, legend="right")

>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2

10.3.2.3 Functional differences

WC from acclimation to FMT2

<<<<<<< HEAD
element_gift %>%
  filter(time_point %in% c("FMT2","Acclimation") & type == "Hot_control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
=======
element_gift %>%
  filter(time_point %in% c("FMT2","Acclimation") & type == "Hot_control") %>%
  select(-time_point, -type) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
  left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>% 
  pivot_longer(-c(Tube_code, time_point),
               names_to = "trait",
               values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
  filter(p_adjust < 0.05) %>%
  remove_rownames() %>% 
  left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>